import numpy as np
import matplotlib.pyplot as plt

# 中文和负号正常显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def analyze_function(f, func_str, x_range):
    """
    全面分析函数的凸性和拐点
    """
    x = np.linspace(x_range[0], x_range[1], 1000)
    y = f(x)
    
    # 使用np.gradient计算导数
    y_prime = np.gradient(y, x)  # 一阶导数
    y_double_prime = np.gradient(y_prime, x)  # 二阶导数（一阶导数的梯度）
    
    # 寻找凸性区间和拐点
    convex_intervals = []  # 下凸区间
    concave_intervals = []  # 上凸区间
    inflection_points = []
    
    current_sign = np.sign(y_double_prime[0])
    start_index = 0
    
    for i in range(1, len(x)):
        sign = np.sign(y_double_prime[i])
        
        if sign != current_sign:
            # 凸性改变，记录拐点
            inflection_points.append((x[i], y[i]))
            
            # 记录前一个区间
            if current_sign > 0:
                convex_intervals.append((x[start_index], x[i]))
            else:
                concave_intervals.append((x[start_index], x[i]))
            
            current_sign = sign
            start_index = i
    
    # 记录最后一个区间
    if current_sign > 0:
        convex_intervals.append((x[start_index], x[-1]))
    else:
        concave_intervals.append((x[start_index], x[-1]))
    
    return x, y, y_prime, y_double_prime, convex_intervals, concave_intervals, inflection_points

# 定义函数
def f(x):
    return 3*x**4 - 4*x**3 + 1

# 执行分析
x, y, y_prime, y_double_prime, convex, concave, inflections = analyze_function(f, "3x⁴ - 4x³ + 1", [-1, 1.5])

# 绘制结果
plt.figure(figsize=(15, 10))

# 1. 原函数
plt.subplot(2, 2, 1)
plt.plot(x, y, 'b-', linewidth=2)
plt.title(r'函数图像: $y = 3x^4 - 4x^3 + 1$')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True, alpha=0.3)

# 2. 一阶导数
plt.subplot(2, 2, 2)
plt.plot(x, y_prime, 'g-', linewidth=2)
plt.title('一阶导数')
plt.xlabel('x')
plt.ylabel("y'")
plt.grid(True, alpha=0.3)

# 3. 二阶导数
plt.subplot(2, 2, 3)
plt.plot(x, y_double_prime, 'r-', linewidth=2)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.title('二阶导数')
plt.xlabel('x')
plt.ylabel("y''")
plt.grid(True, alpha=0.3)

# 4. 凸性分析
plt.subplot(2, 2, 4)
plt.plot(x, y, 'b-', linewidth=2)

# 标记下凸区间
for i, interval in enumerate(convex):
    mask = (x >= interval[0]) & (x <= interval[1])
    if np.any(mask):
        plt.fill_between(x[mask], y[mask], alpha=0.3, color='green', 
                       label='下凸区间' if i == 0 else "")

# 标记上凸区间
for i, interval in enumerate(concave):
    mask = (x >= interval[0]) & (x <= interval[1])
    if np.any(mask):
        plt.fill_between(x[mask], y[mask], alpha=0.3, color='orange', 
                       label='上凸区间' if i == 0 else "")

# 标记拐点
if inflections:
    inflection_x, inflection_y = zip(*inflections)
    plt.plot(inflection_x, inflection_y, 'ro', markersize=8, label='拐点')

plt.title('凸性分析')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("凸性分析结果：")
print("下凸区间 (y'' > 0):", convex)
print("上凸区间 (y'' < 0):", concave)
print("拐点:", inflections)

# 验证结果
print("\n理论验证:")
print("在 x=0 处，y'' =", 36 * 0**2 - 24 * 0, "，函数值 =", f(0))
print("在 x=2/3 处，y'' =", 36*(2/3)**2 - 24*(2/3), "，函数值 =", f(2/3))